home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / polyfitw.pro < prev    next >
Text File  |  1997-07-08  |  2KB  |  103 lines
  1. ; $Id: polyfitw.pro,v 1.2 1996/12/17 20:21:56 ali Exp $
  2.  
  3. FUNCTION POLYFITW,X,Y,W,NDEGREE,YFIT,YBAND,SIGMA,A
  4. ;+
  5. ; NAME:
  6. ;    POLYFITW
  7. ;
  8. ; PURPOSE:
  9. ;    Perform a least-square polynomial fit with optional error estimates.
  10. ;
  11. ; CATEGORY:
  12. ;    Curve fitting.
  13. ;
  14. ; CALLING SEQUENCE:
  15. ;    Result = POLYFITW(X, Y, W, NDegree [, Yfit, Yband, Sigma, A])
  16. ;
  17. ; INPUTS:
  18. ;        X:    The independent variable vector.
  19. ;
  20. ;        Y:    The dependent variable vector.  This vector should be the same 
  21. ;        length as X.
  22. ;
  23. ;        W:    The vector of weights.  This vector should be same length as 
  24. ;        X and Y.
  25. ;
  26. ;     NDegree:    The degree of polynomial to fit.
  27. ;
  28. ; OUTPUTS:
  29. ;    POLYFITW returns a vector of coefficients of length NDegree+1.
  30. ;
  31. ; OPTIONAL OUTPUT PARAMETERS:
  32. ;     Yfit:    The vector of calculated Y's.  Has an error of + or - Yband.
  33. ;
  34. ;    Yband:    Error estimate for each point = 1 sigma.
  35. ;
  36. ;    Sigma:    The standard deviation in Y units.
  37. ;
  38. ;        A:    Correlation matrix of the coefficients.
  39. ;
  40. ; COMMON BLOCKS:
  41. ;    None.
  42. ;
  43. ; SIDE EFFECTS:
  44. ;    None.
  45. ;
  46. ; MODIFICATION HISTORY:
  47. ;    Written by:     George Lawrence, LASP, University of Colorado,
  48. ;            December, 1981.
  49. ;
  50. ;    Adapted to VAX IDL by: David Stern, Jan, 1982.
  51. ;
  52. ;    Weights added, April, 1987,  G. Lawrence
  53. ;
  54. ;-
  55. ON_ERROR,2                      ;RETURN TO CALLER IF AN ERROR OCCURS
  56. N = N_ELEMENTS(X) < N_ELEMENTS(Y) ; SIZE = SMALLER OF X,Y
  57. M = NDEGREE + 1    ; # OF ELEMENTS IN COEFF VEC
  58. ;
  59. A = DBLARR(M,M) ; LEAST SQUARE MATRIX, WEIGHTED MATRIX
  60. B = DBLARR(M)    ; WILL CONTAIN SUM W*Y*X^J
  61. Z = DBLARR(N) + 1.    ; BASIS VECTOR FOR CONSTANT TERM
  62. ;
  63. A[0,0] = TOTAL(W)
  64. B[0] = TOTAL(W*Y)
  65. ;
  66. FOR P = 1,2*NDEGREE DO BEGIN    ; POWER LOOP
  67. Z=Z*X    ; Z IS NOW X^P
  68. IF P LT M THEN B[P] = TOTAL(W*Y*Z)    ; B IS SUM W*Y*X^J
  69. SUM = TOTAL(W*Z)
  70. FOR J = 0 > (P-NDEGREE), NDEGREE < P DO A[J,P-J] = SUM
  71.     END ; END OF P LOOP, CONSTRUCTION OF A AND B
  72. ;
  73. A = INVERT(A)
  74. ;
  75. C = FLOAT(B # A)
  76. ;
  77. IF ( N_PARAMS(0) LE 3) THEN RETURN,C    ; EXIT IF NO ERROR ESTIMATES
  78. ;
  79. ; COMPUTE OPTIONAL OUTPUT PARAMETERS.
  80. ;
  81. YFIT = FLTARR(N)+C[0]    ; ONE-SIGMA ERROR ESTIMATES, INIT
  82. FOR K = 1, NDEGREE DO YFIT = YFIT + C[K]*(X^K)    ; SUM BASIS VECTORS
  83. ;
  84. VAR = TOTAL((YFIT-Y)^2 )/(N-M)    ; VARIANCE ESTIMATE, UNBIASED
  85. ;
  86. ;
  87. SIGMA=SQRT(VAR)
  88. YBAND = FLTARR(N) + A[0,0]
  89. Z=FLTARR(N)+1.
  90. FOR P=1,2*NDEGREE DO BEGIN    ; COMPUTE CORRELATED ERROR ESTIMATES ON Y
  91.     Z = Z*X        ; Z IS NOW X^P
  92.     SUM = 0.
  93.     FOR J=0 > (P - NDEGREE), NDEGREE  < P DO SUM = SUM + A[J,P-J]
  94.     YBAND = YBAND + SUM * Z    ; ADD IN ALL THE ERROR SOURCES
  95. ENDFOR    ; END OF P LOOP
  96.     YBAND = YBAND*VAR
  97.     YBAND = SQRT( YBAND )
  98. RETURN,C
  99. END
  100.  
  101.  
  102.  
  103.